home *** CD-ROM | disk | FTP | other *** search
/ IRIX 6.2 Development Libraries / SGI IRIX 6.2 Development Libraries.iso / dist / complib.idb / usr / share / catman / p_man / cat3 / complib / ctrsen.z / ctrsen
Text File  |  1996-03-14  |  9KB  |  265 lines
  1.  
  2.  
  3.  
  4. CCCCTTTTRRRRSSSSEEEENNNN((((3333FFFF))))                                                          CCCCTTTTRRRRSSSSEEEENNNN((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      CTRSEN - reorder the Schur factorization of a complex matrix A =
  10.      Q*T*Q**H, so that a selected cluster of eigenvalues appears in the
  11.      leading positions on the diagonal of the upper triangular matrix T, and
  12.      the leading columns of Q form an orthonormal basis of the corresponding
  13.      right invariant subspace
  14.  
  15. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  16.      SUBROUTINE CTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP,
  17.                         WORK, LWORK, INFO )
  18.  
  19.          CHARACTER      COMPQ, JOB
  20.  
  21.          INTEGER        INFO, LDQ, LDT, LWORK, M, N
  22.  
  23.          REAL           S, SEP
  24.  
  25.          LOGICAL        SELECT( * )
  26.  
  27.          COMPLEX        Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
  28.  
  29. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  30.      CTRSEN reorders the Schur factorization of a complex matrix A = Q*T*Q**H,
  31.      so that a selected cluster of eigenvalues appears in the leading
  32.      positions on the diagonal of the upper triangular matrix T, and the
  33.      leading columns of Q form an orthonormal basis of the corresponding right
  34.      invariant subspace.
  35.  
  36.      Optionally the routine computes the reciprocal condition numbers of the
  37.      cluster of eigenvalues and/or the invariant subspace.
  38.  
  39.  
  40. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  41.      JOB     (input) CHARACTER*1
  42.              Specifies whether condition numbers are required for the cluster
  43.              of eigenvalues (S) or the invariant subspace (SEP):
  44.              = 'N': none;
  45.              = 'E': for eigenvalues only (S);
  46.              = 'V': for invariant subspace only (SEP);
  47.              = 'B': for both eigenvalues and invariant subspace (S and SEP).
  48.  
  49.      COMPQ   (input) CHARACTER*1
  50.              = 'V': update the matrix Q of Schur vectors;
  51.              = 'N': do not update Q.
  52.  
  53.      SELECT  (input) LOGICAL array, dimension (N)
  54.              SELECT specifies the eigenvalues in the selected cluster. To
  55.              select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
  56.  
  57.      N       (input) INTEGER
  58.              The order of the matrix T. N >= 0.
  59.  
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. CCCCTTTTRRRRSSSSEEEENNNN((((3333FFFF))))                                                          CCCCTTTTRRRRSSSSEEEENNNN((((3333FFFF))))
  71.  
  72.  
  73.  
  74.      T       (input/output) COMPLEX array, dimension (LDT,N)
  75.              On entry, the upper triangular matrix T.  On exit, T is
  76.              overwritten by the reordered matrix T, with the selected
  77.              eigenvalues as the leading diagonal elements.
  78.  
  79.      LDT     (input) INTEGER
  80.              The leading dimension of the array T. LDT >= max(1,N).
  81.  
  82.      Q       (input/output) COMPLEX array, dimension (LDQ,N)
  83.              On entry, if COMPQ = 'V', the matrix Q of Schur vectors.  On
  84.              exit, if COMPQ = 'V', Q has been postmultiplied by the unitary
  85.              transformation matrix which reorders T; the leading M columns of
  86.              Q form an orthonormal basis for the specified invariant subspace.
  87.              If COMPQ = 'N', Q is not referenced.
  88.  
  89.      LDQ     (input) INTEGER
  90.              The leading dimension of the array Q.  LDQ >= 1; and if COMPQ =
  91.              'V', LDQ >= N.
  92.  
  93.      W       (output) COMPLEX
  94.              The reordered eigenvalues of T, in the same order as they appear
  95.              on the diagonal of T.
  96.  
  97.      M       (output) INTEGER
  98.              The dimension of the specified invariant subspace.  0 <= M <= N.
  99.  
  100.      S       (output) REAL
  101.              If JOB = 'E' or 'B', S is a lower bound on the reciprocal
  102.              condition number for the selected cluster of eigenvalues.  S
  103.              cannot underestimate the true reciprocal condition number by more
  104.              than a factor of sqrt(N). If M = 0 or N, S = 1.  If JOB = 'N' or
  105.              'V', S is not referenced.
  106.  
  107.      SEP     (output) REAL
  108.              If JOB = 'V' or 'B', SEP is the estimated reciprocal condition
  109.              number of the specified invariant subspace. If M = 0 or N, SEP =
  110.              norm(T).  If JOB = 'N' or 'E', SEP is not referenced.
  111.  
  112.      WORK    (workspace) COMPLEX array, dimension (LWORK)
  113.              If JOB = 'N', WORK is not referenced.
  114.  
  115.      LWORK   (input) INTEGER
  116.              The dimension of the array WORK.  If JOB = 'N', LWORK >= 1; if
  117.              JOB = 'E', LWORK = M*(N-M); if JOB = 'V' or 'B', LWORK >=
  118.              2*M*(N-M).
  119.  
  120.      INFO    (output) INTEGER
  121.              = 0:  successful exit
  122.              < 0:  if INFO = -i, the i-th argument had an illegal value
  123.  
  124.  
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. CCCCTTTTRRRRSSSSEEEENNNN((((3333FFFF))))                                                          CCCCTTTTRRRRSSSSEEEENNNN((((3333FFFF))))
  137.  
  138.  
  139.  
  140. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  141.      CTRSEN first collects the selected eigenvalues by computing a unitary
  142.      transformation Z to move them to the top left corner of T. In other
  143.      words, the selected eigenvalues are the eigenvalues of T11 in:
  144.  
  145.                    Z'*T*Z = ( T11 T12 ) n1
  146.                             (  0  T22 ) n2
  147.                                n1  n2
  148.  
  149.      where N = n1+n2 and Z' means the conjugate transpose of Z. The first n1
  150.      columns of Z span the specified invariant subspace of T.
  151.  
  152.      If T has been obtained from the Schur factorization of a matrix A =
  153.      Q*T*Q', then the reordered Schur factorization of A is given by A =
  154.      (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the
  155.      corresponding invariant subspace of A.
  156.  
  157.      The reciprocal condition number of the average of the eigenvalues of T11
  158.      may be returned in S. S lies between 0 (very badly conditioned) and 1
  159.      (very well conditioned). It is computed as follows. First we compute R so
  160.      that
  161.  
  162.                             P = ( I  R ) n1
  163.                                 ( 0  0 ) n2
  164.                                   n1 n2
  165.  
  166.      is the projector on the invariant subspace associated with T11.  R is the
  167.      solution of the Sylvester equation:
  168.  
  169.                            T11*R - R*T22 = T12.
  170.  
  171.      Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote the
  172.      two-norm of M. Then S is computed as the lower bound
  173.  
  174.                          (1 + F-norm(R)**2)**(-1/2)
  175.  
  176.      on the reciprocal of 2-norm(P), the true reciprocal condition number.  S
  177.      cannot underestimate 1 / 2-norm(P) by more than a factor of sqrt(N).
  178.  
  179.      An approximate error bound for the computed average of the eigenvalues of
  180.      T11 is
  181.  
  182.                             EPS * norm(T) / S
  183.  
  184.      where EPS is the machine precision.
  185.  
  186.      The reciprocal condition number of the right invariant subspace spanned
  187.      by the first n1 columns of Z (or of Q*Z) is returned in SEP.  SEP is
  188.      defined as the separation of T11 and T22:
  189.  
  190.                         sep( T11, T22 ) = sigma-min( C )
  191.  
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. CCCCTTTTRRRRSSSSEEEENNNN((((3333FFFF))))                                                          CCCCTTTTRRRRSSSSEEEENNNN((((3333FFFF))))
  203.  
  204.  
  205.  
  206.      where sigma-min(C) is the smallest singular value of the
  207.      n1*n2-by-n1*n2 matrix
  208.  
  209.         C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
  210.  
  211.      I(m) is an m by m identity matrix, and kprod denotes the Kronecker
  212.      product. We estimate sigma-min(C) by the reciprocal of an estimate of the
  213.      1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) cannot
  214.      differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
  215.  
  216.      When SEP is small, small changes in T can cause large changes in the
  217.      invariant subspace. An approximate bound on the maximum angular error in
  218.      the computed right invariant subspace is
  219.  
  220.                          EPS * norm(T) / SEP
  221.  
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.